//this code deals with meshing isosurfaces to procedurally generate
//smooth surfaces using the surface nets method

//sphere isosurface
//an implicit surface represented by a field function
var sphere_isosurface = function() {
    this.center = new THREE.Vector3(0,0,0);
    this.radius = 5;
    var v1 = new THREE.Vector3().copy(this.center);
    var v2 = new THREE.Vector3().copy(this.center);
    this.aabb = new THREE.Box3(v1.addScalar(-this.radius),v2.addScalar(this.radius));
    this.transform = new THREE.Matrix4();

    var t = this;

    //returns the normal at some point
    this.normal = function(x,y,z) {
        
    }

    //given a point it returns an array of densities with size of num_points
    //TODO: transform the given x,y,z using the inverse of our
    //transformation matrix
    this.density = function(x, y, z, dz, num_points) {
        num_points = num_points || 10;
        var densities = []
        for(var i=0; i<num_points; ++i) {
            var sq_dist = x*x + y*y + z*z;
            var sq_rad = t.radius*t.radius;
            densities[i] = sq_dist - sq_rad;
            z += dz;
        }
        return densities;
    }
}

/**
 * SurfaceNets in JavaScript
 *
 * Written by Mikola Lysenko (C) 2012
 *
 * MIT License
 *
 * Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
 */
"use strict";

//modified from original source:
//https://github.com/mikolalysenko/isosurface/blob/master/lib/surfacenets.js


//Precompute edge table, like Paul Bourke does.
// This saves a bit of time when computing the centroid of each boundary cell
var cube_edges = new Int32Array(24)
  , edge_table = new Int32Array(256);
(function() {

  //Initialize the cube_edges table
  // This is just the vertex number of each cube
  var k = 0;
  for(var i=0; i<8; ++i) {
    for(var j=1; j<=4; j<<=1) {
      var p = i^j;
      if(i <= p) {
        cube_edges[k++] = i;
        cube_edges[k++] = p;
      }
    }
  }

  //Initialize the intersection table.
  //  This is a 2^(cube configuration) ->  2^(edge configuration) map
  //  There is one entry for each possible cube configuration, and the output is a 12-bit vector enumerating all edges crossing the 0-level.
  for(var i=0; i<256; ++i) {
    var em = 0;
    for(var j=0; j<24; j+=2) {
      var a = !!(i & (1<<cube_edges[j]))
        , b = !!(i & (1<<cube_edges[j+1]));
      em |= a !== b ? (1 << (j >> 1)) : 0;
    }
    edge_table[i] = em;
  }
})();

//Internal buffer, this may get resized at run time
var buffer = new Array(4096);
(function() {
  for(var i=0; i<buffer.length; ++i) {
    buffer[i] = 0;
  }
})();


//dims: resolution of the isosurface
//potential: scalar valued potential function taking 3 coords
//bounds: a pair of 3d vectors [lo, hi] giving bounds of potential to sample
var surfaceNets = function(dims, potential, bounds) {
    if(!bounds) {
        bounds = [[0,0,0],dims];
    }
    
    var scale     = [0,0,0];
    var shift     = [0,0,0];
    for(var i=0; i<3; ++i) {
        scale[i] = (bounds[1][i] - bounds[0][i]) / dims[i];
        shift[i] = bounds[0][i];
    }
    
    var vertices = []
    , faces = []
    , n = 0
    , x = [0, 0, 0]
    , R = [1, (dims[0]+1), (dims[0]+1)*(dims[1]+1)]
    , grid = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    , buf_no = 1;
    
    
    //Resize buffer if necessary 
    if(R[2] * 2 > buffer.length) {
        var ol = buffer.length;
        buffer.length = R[2] * 2;
        while(ol < buffer.length) {
            buffer[ol++] = 0;
        }
    }
    
    //March over the voxel grid
    for(x[2]=0; x[2]<dims[2]-1; ++x[2], n+=dims[0], buf_no ^= 1, R[2]=-R[2]) {
        
        //m is the pointer into the buffer we are going to use.  
        //This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
        //The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
        var m = 1 + (dims[0]+1) * (1 + buf_no * (dims[1]+1));
        
        for(x[1]=0; x[1]<dims[1]-1; ++x[1], ++n, m+=2)
            for(x[0]=0; x[0]<dims[0]-1; ++x[0], ++n, ++m) {
                
                //Read in 8 field values around this vertex and store them in an array
                //Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
                var mask = 0, g = 0;
                for(var k=0; k<2; ++k)
                    for(var j=0; j<2; ++j)      
                        for(var i=0; i<2; ++i, ++g) {
                            var p = potential(
                                scale[0]*(x[0]+i)+shift[0],
                                scale[1]*(x[1]+j)+shift[1],
                                scale[2]*(x[2]+k)+shift[2]);
                            grid[g] = p;
                            mask |= (p < 0) ? (1<<g) : 0;
                        }
                
                //Check for early termination if cell does not intersect boundary
                if(mask === 0 || mask === 0xff) {
                    continue;
                }
                
                //Sum up edge intersections
                var edge_mask = edge_table[mask]
                , v = [0.0,0.0,0.0]
                , e_count = 0;
                
                //For every edge of the cube...
                for(var i=0; i<12; ++i) {
                    
                    //Use edge mask to check if it is crossed
                    if(!(edge_mask & (1<<i))) {
                        continue;
                    }
                    
                    //If it did, increment number of edge crossings
                    ++e_count;
                    
                    //Now find the point of intersection
                    var e0 = cube_edges[ i<<1 ]       //Unpack vertices
                    , e1 = cube_edges[(i<<1)+1]
                    , g0 = grid[e0]                 //Unpack grid values
                    , g1 = grid[e1]
                    , t  = g0 - g1;                 //Compute point of intersection
                    if(Math.abs(t) > 1e-6) {
                        t = g0 / t;
                    } else {
                        continue;
                    }
                    
                    //Interpolate vertices and add up intersections (this can be done without multiplying)
                    for(var j=0, k=1; j<3; ++j, k<<=1) {
                        var a = e0 & k
                        , b = e1 & k;
                        if(a !== b) {
                            v[j] += a ? 1.0 - t : t;
                        } else {
                            v[j] += a ? 1.0 : 0;
                        }
                    }
                }
                
                //Now we just average the edge intersections and add them to coordinate
                var s = 1.0 / e_count;
                for(var i=0; i<3; ++i) {
                    v[i] = scale[i] * (x[i] + s * v[i]) + shift[i];
                }
                
                //Add vertex to buffer, store pointer to vertex index in buffer
                buffer[m] = vertices.length;
                vertices.push(v);
                
                //Now we need to add faces together, to do this we just loop over 3 basis components
                for(var i=0; i<3; ++i) {
                    //The first three entries of the edge_mask count the crossings along the edge
                    if(!(edge_mask & (1<<i)) ) {
                        continue;
                    }
                    
                    // i = axes we are point along.  iu, iv = orthogonal axes
                    var iu = (i+1)%3
                    , iv = (i+2)%3;
                    
                    //If we are on a boundary, skip it
                    if(x[iu] === 0 || x[iv] === 0) {
                        continue;
                    }
                    
                    //Otherwise, look up adjacent edges in buffer
                    var du = R[iu]
                    , dv = R[iv];
                    
                    //Remember to flip orientation depending on the sign of the corner.
                    if(mask & 1) {
                        faces.push([buffer[m],    buffer[m-du],    buffer[m-dv]]);
                        faces.push([buffer[m-dv], buffer[m-du],    buffer[m-du-dv]]);
                    } else {
                        faces.push([buffer[m],    buffer[m-dv],    buffer[m-du]]);
                        faces.push([buffer[m-du], buffer[m-dv],    buffer[m-du-dv]]);
                    }
                }
            }
    }
    
    //All done!  Return the result
    //return { positions: vertices, cells: faces };
    return { vertices: vertices, faces: faces };
};

//takes in an array containing the data from a field function along
//with the resolution
var surfaceNets2 = function(data, dims) {
    var vertices = []
    , faces = []
    , n = 0
    , x = new Int32Array(3)
    , R = new Int32Array([1, (dims[0]+1), (dims[0]+1)*(dims[1]+1)])
    , grid = new Float32Array(8)
    , buf_no = 1;
    
    //Resize buffer if necessary 
    if(R[2] * 2 > buffer.length) {
        buffer = new Int32Array(R[2] * 2);
    }
    
    //March over the voxel grid
    for(x[2]=0; x[2]<dims[2]-1; ++x[2], n+=dims[0], buf_no ^= 1, R[2]=-R[2]) {
        //m is the pointer into the buffer we are going to use.  
        //This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
        //The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
        var m = 1 + (dims[0]+1) * (1 + buf_no * (dims[1]+1));
        
        for(x[1]=0; x[1]<dims[1]-1; ++x[1], ++n, m+=2)
            for(x[0]=0; x[0]<dims[0]-1; ++x[0], ++n, ++m) {
                
                //Read in 8 field values around this vertex and store them in an array
                //Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
                var mask = 0, g = 0, idx = n;
                for(var k=0; k<2; ++k, idx += dims[0]*(dims[1]-2))
                    for(var j=0; j<2; ++j, idx += dims[0]-2)      
                        for(var i=0; i<2; ++i, ++g, ++idx) {
                            var p = data[idx];
                            grid[g] = p;
                            mask |= (p < 0) ? (1<<g) : 0;
                        }
                
                //Check for early termination if cell does not intersect boundary
                if(mask === 0 || mask === 0xff) {
                    continue;
                }
                
                //Sum up edge intersections
                var edge_mask = edge_table[mask]
                , v = [0.0,0.0,0.0]
                , e_count = 0;
                
                //For every edge of the cube...
                for(var i=0; i<12; ++i) {
                    
                    //Use edge mask to check if it is crossed
                    if(!(edge_mask & (1<<i))) {
                        continue;
                    }
                    
                    //If it did, increment number of edge crossings
                    ++e_count;
                    
                    //Now find the point of intersection
                    var e0 = cube_edges[ i<<1 ]       //Unpack vertices
                    , e1 = cube_edges[(i<<1)+1]
                    , g0 = grid[e0]                 //Unpack grid values
                    , g1 = grid[e1]
                    , t  = g0 - g1;                 //Compute point of intersection
                    if(Math.abs(t) > 1e-6) {
                        t = g0 / t;
                    } else {
                        continue;
                    }
                    
                    //Interpolate vertices and add up intersections (this can be done without multiplying)
                    for(var j=0, k=1; j<3; ++j, k<<=1) {
                        var a = e0 & k
                        , b = e1 & k;
                        if(a !== b) {
                            v[j] += a ? 1.0 - t : t;
                        } else {
                            v[j] += a ? 1.0 : 0;
                        }
                    }
                }
                
                //Now we just average the edge intersections and add them to coordinate
                var s = 1.0 / e_count;
                for(var i=0; i<3; ++i) {
                    v[i] = x[i] + s * v[i];
                }
                
                //Add vertex to buffer, store pointer to vertex index in buffer
                buffer[m] = vertices.length;
                vertices.push(v);
                
                //Now we need to add faces together, to do this we just loop over 3 basis components
                for(var i=0; i<3; ++i) {
                    //The first three entries of the edge_mask count the crossings along the edge
                    if(!(edge_mask & (1<<i)) ) {
                        continue;
                    }
                    
                    // i = axes we are point along.  iu, iv = orthogonal axes
                    var iu = (i+1)%3
                    , iv = (i+2)%3;
                    
                    //If we are on a boundary, skip it
                    if(x[iu] === 0 || x[iv] === 0) {
                        continue;
                    }
                    
                    //Otherwise, look up adjacent edges in buffer
                    var du = R[iu]
                    , dv = R[iv];
                    
                    //Remember to flip orientation depending on the sign of the corner.
                    if(mask & 1) {
                        faces.push([buffer[m], buffer[m-du], buffer[m-du-dv], buffer[m-dv]]);
                    } else {
                        faces.push([buffer[m], buffer[m-dv], buffer[m-du-dv], buffer[m-du]]);
                    }
                }
            }
    }
    
    //All done!  Return the result
    return { vertices: vertices, faces: faces };
};


function volume_cache(f) {
    var cached = null;
    return function() {
        if(cached === null) { 
            cached = f();
        }
        return cached;
    }
}

function make_volume(dims, f) {
    return volume_cache(function() {
        var res = new Array(3);
        for(var i=0; i<3; ++i) {
            res[i] = 2 + Math.ceil((dims[i][1] - dims[i][0]) / dims[i][2]);
        }
        var volume = new Float32Array(res[0] * res[1] * res[2])
        , n = 0;
        for(var k=0, z=dims[2][0]-dims[2][2]; k<res[2]; ++k, z+=dims[2][2])
            for(var j=0, y=dims[1][0]-dims[1][2]; j<res[1]; ++j, y+=dims[1][2])
                for(var i=0, x=dims[0][0]-dims[0][2]; i<res[0]; ++i, x+=dims[0][2], ++n) {
                    volume[n] = f(x,y,z);
                }
        return {data: volume, dims:res};
    });
}

//an isosurface that is polygonized using the surface nets method
var isosurface = function(options) {
    var t = this;

    //options
    this.color = new THREE.Color(0x00edff);
    this.field = null; //field function
    this.dims = null; //dimensions
    this.volume = null; //actual volume data and dimensions passed ro surfaceNets2
    this.potential = null; //potential function
    this.cache = {};
    this.mesh = null;
    this.wiremesh = null;
    this.geo = null;
    this.mat = new THREE.MeshNormalMaterial(); 
    this.wiremat = new THREE.MeshBasicMaterial({color: this.color});
    this.object = null;

    this.set_options = function(opts) {
        for (o in opts) {
            var val = opts[o];
            if(val === undefined) {
                console.warn('isosurface: undefined value ' + val + ' passed as option.');
                continue;
            }

            if(o in t) {
                var curVal = t[o];
                if(curVal instanceof THREE.Color) {
                    curVal.set(val);
                } else if(curVal instanceof THREE.Vector3 && val instanceof THREE.Vector3) {
                    curVal.copy(val);
                } else {
                    t[o] = val;
                }
            }
        }
    };

    this.set_options(options);

    this.build_volume = function() {
        if(!t.field || !t.dims) {
            console.warn('isosurface: cannot build volume without field function and dimensions.');
        }

        t.volume = make_volume(t.dims, t.field)();
    }

    this.build_mesh = function() {
        var geo;
        if(!t.potential && !t.volume) {
            console.warn('isosurface: cannot build mesh without volume data or potential function');
        }

        var sn = null;
        if(t.potential) {
            sn = surfaceNets([64,64,64],
                             t.potential
                             ,[[-10,-10,-10],[10,10,10]]);  
        } else if(t.volume) {
            sn = surfaceNets2(t.volume.data, t.volume.dims);
        }

        geo = new THREE.Geometry();
        for(var i=0; i<sn.vertices.length; ++i) {
            var v = sn.vertices[i];
            geo.vertices.push(new THREE.Vector3(v[0], v[1], v[2]));
        }
        
        for(var i=0; i<sn.faces.length; i++) {
            var f = sn.faces[i];
            if(f.length === 3) {
                geo.faces.push(new THREE.Face3(f[0], f[1], f[2]));
            }
            else if(f.length === 4) {
                geo.faces.push(new THREE.Face4(f[0], f[1], f[2], f[3]));
            }
        }

        geo.computeFaceNormals();
        geo.verticesNeedUpdate = true;
        geo.elementsNeedUpdate = true;
        geo.normalsNeedUpdate = true;

        geo.computeBoundingBox();
        geo.computeBoundingSphere();
        
        t.geo = geo;
        t.mesh = new THREE.Mesh(t.geo, t.mat);
        t.wiremesh = new THREE.Mesh(t.geo, t.wiremat);

        t.object = new THREE.Object3D();
        t.object.add(t.mesh);
        t.object.add(t.wiremesh);
        return t.object;
    }

    //check to see if a point it intersecting this isosurface
    //for now I'll just raycast against the geometry surfaceNets makes
    this.intersects = function (point) {
        if(!t.volume) {
            console.warn('isosurface: cannot calculate intersection without volume data.');
        }
        //do tranforms to the point based on object3d tranform
        //information and test point against field function to get a
        //density value
    }

}

//TODO: look into FRep and R-Functions
